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Abstract 

Brownian dynamics simulations were carried out to study wave spectra of two-dimensional dusty 
plasma liquids and solids for a wide range of wavelengths. The existence of a longitudinal dust 
thermal mode was confirmed in simulations, and a cutoff wavenumber in the transverse mode 
was measured. Dispersion relations, resulting from simulations, were compared with those from 
analytical theories, such as the random-phase approximation (RPA), quasi- localized charged ap- 
proximation (QLCA), and harmonic approximation (HA). An overall good agreement between the 
QLCA and simulations was found for wide ranges of states and wavelengths after taking into ac- 
count the direct thermal effect in the QLCA, while for the RPA and HA good agreement with 
simulations were found in the high and low temperature limits, respectively. 

PACS numbers: 52.25.Fi, 52.27.Gr, 52.27.Lw 
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I. INTRODUCTION 



A laboratory-generated dusty plasma is a suspension of micron-sized particles immersed 
in a usual plasma with ions, electrons and neutral gas molecules Dust particles 

acquire a few thousand of electron charges by absorbing the surrounding electrons and 
ions, and they consequently interact with each other via a dynamically-screened Coulomb 
potential [4j, while undergoing Brownian motions due to frequent collisions, mainly with the 
neutral gas molecules. When the interaction potential energy between charged dust particles 
significantly exceeds their kinetic energy, they become strongly coupled and, consequently, 
they can form ordered structures characteristic of a liquid or solid state. Such structures 
are commonly referred to as strongly coupled dusty plasmas (SCDPs). 

Two-dimensional (2D) SCDPs have become particularly favored in recent laboratory ex- 
periments, because in 2D geometry the complication due to ion wake effect (see, e. g., [4j and 
references therein) should be absent and the particle interaction can be well approximated 



by the Yukawa potential 



BB. 



Of particular interest in 2D dusty plasmas are their collec- 



tive and dynamical properties, such as longitudinal and transverse wave modes, which have 
)een studied extensively over the past decade in experiments 3, S S H, [ill Q] , theories 
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14. 



15 



16 



17 



18 



19j , and numerical simulations 



18 



20 



2l| . On the experimental side, 



externally^excited (longitudinal) dust lattice waves (DLWs) were first observed by Homann 
et al. 



I. 0, 



|8J, who found their dispersion relation to be in good agreement with the theoreti- 
cal prediction of Melands0 14j. Next, thermally excited phonon spectrum in a 2D plasma 
crystal was observed in an experiment by Nunomura et al. |9], who demonstrated a good 



agreement of this spectrum with theory 



15 



16| for both the longitudinal and transverse 



modes in the entire first Brillouin zone. That experiment was later extended by Zhdanov 
et al. [h]], who studied the polarization of wave modes in a 2D plasma crystal, and by 
Nunomura et al. llj], who studied the wave spectra in both liquid and solid states for a 
wide range of wavelengths. More recently, Nosenko et al. [12] have measured experimentally 
a critical cutoff wave number for shear waves in 2D dusty plasma liquids. 

On the theoretical side, besides the above mentioned DLW theories of Peeters and Wu 
Q, Melandwfl, Dub^, and Wang et al. Q, a s wen as the standard dnst aeon Stl e 
wave (DAW) theory {22, 23] along with its 2D derivatives 24, 25], there have been many 
other studies of SCDPs. For example, a semi-analytic approximation was used to study 



2 



wave propagation in 2D SCDPs [15j, |25J. The quasi-localized charge approximation (QLCA) 



261 ] was used to study collective modes and dynamics of both 3D 



221 



and 2D 



dusty plasma liquids. Furthermore, a generalized hydrodynamic (GH) model was adopted 
by Kaw and Sen 2^| to study wave dispersion in 3D dusty plasma liquids, while Murillo and 
coworkers studied collective modes by using both kinetic |3o| and the GH 31, 32] methods. 
The latter method was used, in particular, to study critical wavenumbers for transverse 



modes in a 3D dusty plasma in liquid phase 



20 



21 



33 



3jj. On the other hand, computer simulations 



34|, as an essential supplement to real experiments, have played important 



roles in validating various analytical theories and in explaining experimental observations. 
Both Molecular dynamics (MD) and the particle-in-cell (PIC) simulations were carried out 



by Winske et al. 



33] to study longitudinal wave dispersion in one-dimensional systems. 



271 ] including the 



Results were compared with the DAW mode and with the 3D QLCA 
strong-coupling effects, and an agreement with the QLCA was found only at very long 
wavelengths due to the small system size. Later on, Ohta and Hamaguchi 34j studied both 
longitudinal and transverse dispersion relations in 3D dust y p lasma liquids, and results were 



compared with both the QLCA 27, 



29 



311 ] . It was found that, for the 



and the GH model 

longitudinal mode, both theories were in close agreement with simulations in the entire first 
Brillouin zone, whereas, for the transverse mode, good agreements with simulations were 
found for both theories, except in the very long wavelength region, where the GH model 
successfully predicts a cutoff but the QLCA does not. This deficiency of the QLCA can be 
rectified by introducing a phenomenological damping which accounts for the diffusional and 
other damping effects, as was shown more recently in the work by Kalman et al. [18l| . where 
the QLCA was extended to study 2D dusty plasma liquids, and the results were critically 
compared with MD simulations. Good agreements with simulations were found for both the 
longitudinal and transverse modes, mostly in the first Brillouin zone. (Please see Ref. 35] 
for a nice review of recent development about collective dynamics in 2D and 3D Yukawa 
liquids, based mainly on QLCA and computer simulations.) 

However, it should be noted that, first, most of the above simulations 

considered dusty plasmas in the liquid state, and the obtained dispersion relations 



18 



20 



21 



33, 



were limited to the first Brillouin zone. To the best of our knowledge, there are no simu- 
lations verifying the above analytical theories for a wider range of wavelengths, as well as 
for dusty plasmas in non-ideal gaseous, or in solid states. Secondly, no simulations have 
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been conducted to show a transition from the random phase approximation (RPA) based 
DAW to the harmonic approximation (HA) based DLW when dusty plasma goes from a 
high temperature liquid (or non-ideal gaseous) state to a low temperature crystalline state. 
Third, an analytical theory based on the RPA {23] predicted a so-called dust thermal wave 
(DTW) due to the direct thermal effect of dust, and its existence was verified in a subsequent 



experiment 



ll). However, the original theory 23] applies only to weakly coupled systems, 



and it would be interesting to study, both analytically and via simulation, how the DTW 
behaves in strongly coupled systems. 

To fulfill these goals, we perform here computer simulations and study the resulting wave 
spectra of 2D dusty plasmas for a wide range of wavelengths and system states. Results are 
compared with the available analytical theories for the 2D geometry, such as, RPA, QLCA 
and HA. In addition, two extensions of the standard QLCA are discussed, one including the 
direct thermal effect, and the other implementing a critical wave cutoff in the transverse 



18 



351 ] . The role of damping effects on 



mode, based on the method of Kalman et al. 
collective modes is also investigated in the simulation. The remaining part of the paper is 
organized as follows. Details of our simulations are presented in Sec. II, which is followed 
by a presentation and discussion of the results in Sec. Ill, including brief derivations and 
reviews of various analytical theories. Concluding remarks are given in Sec. IV. 



II. SIMULATION 
A. Algorithm 



Our simulation is based on the Brownian dynamics (BD) method [36l . 1371 . l38j . which 
may be regarded as a generalization of the standard MD method. Namely, while the MD 
simulation is based on Newton's equations of motion, the BD method is based on their 
generalization in the form of Langevin equation (and its integral), viz., 
d 

— r = v 

dt 

= -uv+-F + A(t), (1) 
dt m 

where, as usual, m, v and r are, respectively, the mass, velocity and position of a Brow- 
nian particle, and F is the systematic (deterministic) force coming from the inter-particle 
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interactions within the system and, possibly, from external force fields. What is different 
from Newton's equations is the appearance of dynamic friction, — vv, and the Brownian ac- 
celeration, A(£), which represent complementing effects of a single, sub-scale phenomenon: 
numerous, frequent collisions of Brownian particle with molecules in the medium. While 
the former represents the average effect of these collisions, the latter represents fluctuations 
due to discreteness of the collisions, and is generally assumed to be well represented by a 
delta-correlated Gaussian white noise. They are both related to the medium temperature 
through a fluctuation-dissipation theorem. 

The Langevin equations Eq. (CD) may be numerically integrated in a manner similar to the 
integration of Newton's equations in the MD method, and such a technique is generally called 
Brownian dynamics. (Note that several different names are used in literature to designate 
methods for numerical integration of Langevin equation, for example, Brownian Dynamics, 
Langevin Dynamics, or Langevin Molecular Dynamics, depending on the background, area, 
or preferences of different researchers. Those names may refer to the same technique that we 
discuss here, or they may involve subtle differences in derivations of the simulation formulae, 
and/or in implementations of simulation. We follow here definitions given by Allen and 



Tildesley 



36[.) The advantages of BD here are as follows. Firstly it brings the simulation 



closer to real dusty plasma experiments by taking into account both the Brownian motion 
and damping effect self-consistently. And secondly it simplifies the simulation in a sense 
that no external thermostat is need. 

for BD 



heat conduction 



38 



We employ here the 5th order Gear-like predictor-corrector algorithm 37, 
simulation, which was used successfully in simulating shock wave propagation [39, 

43], and diffusion processes in SCDPs [44J]. When compared with several 



40 



41 



3, 



other popular methods for BD simulation, such as Euler-like, Beeman-like and Verlet-like 
methods 36|, the present method can cover a wider range of friction coefficients u, and is 



particularly reliable in the low-damping regime while exhibiting higher-order accuracy, good 
stability, and negligible drift on long time scales. Therefore, the present algorithm appears 
to be especially suitable to simulate dusty plasmas, which are often slightly damped. 
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FIG. 1: (Color online) Radial distribution function g{r) for n = 1 and different V. 



B. Calculation of wave spectra 



Typically, iV = 4000 dust particles are simulated in a square with periodic boundary 
conditions. Particles interact with each other via pairwise Yukawa potentials 6j of the 
form <j){r) = (e 2 Z%/r) exp (— r/\r>), where e is the elementary charge, is the average 
number of charges on each dust particle, r the inter-particle distance, and \p the Debye 
screening length. Note here that the Yukawa interaction is not merely an assumption, but 
is rather consistent with theoretical derivations, described in the following section. Such a 
system can be fully characterized by three parameters 45|]: the Coulomb coupling parameter, 
T = (eZd) 2 /(aT), the screening parameter, n = a/Xo, and the damping rate, v/uj p d, due 
to the neutral gas, where T is the system temperature (in energy units), a = (nado) 1 ^ 2 
the Wigner-Seitz radius, and a^o the equilibrium dust number density. The dust plasma 
frequency is defined by u p d = [2(eZ d ) 2 / (ma 3 )] 1 ' 2 , where m is the dust particle mass. We 
shall use k = 1 throughout our simulation and in subsequent discussions, because k ~ 1 is 
the most typical value found in experiments and, moreover, no substantial effects of varying 
Kj are expected in the features to be discussed in the following. 

Initially, dust particles are randomly placed in the square. The system comes to an 
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equilibrium after a period which is proportional to 1/u, and is usually about 10/u. The 
time step used in simulations is 5t = 0.02^"^ for T > 10 and 5t = O.Ola;^ 1 for T < 10. After 
the system reaches an equilibrium, radial distribution function g{r) and wave spectra are 
calculated. 

Figure [T] shows examples of g(r) for k = 1, with different T values. These results will be 
used as an input below in evaluating the QLCA dispersion relations. 

The longitudinal, C(k, lu), and transverse, T(k, lu), wave spectra are determined by means 
of the current- current correlation functions in the longitudinal and transverse directions, 



respectively. We follow Ref. [461 ] . in which the current density of the system is defined by 

N 



jMH^Ev^r-r^)], (2) 

* i=i 

with the Fourier transform of its cartesian component a given by 
1 N 

j a ( k ,*) = — ^ a (tK k - r <w (3) 

i=l 

where rj(t) and Vj(t) are, respectively, position and velocity of the ith particle at time t, 
with a being x or y. Assuming that waves propagate along the x direction, i.e., k = {k, 0}, 
one defines 

J L (k,t) = (j*(k,t)j x (k,0)), 

j T (k,t) = (j;(k,t) h (k,o)), 

to be, respectively, longitudinal and transverse current auto-correlations. Here, the asterisk 
designates complex conjugation, while the angular brackets indicate an ensemble average. 
The longitudinal and transverse wave spectra are then obtained by Fourier transforms of 



the corresponding current auto-correlations 46] 

/+oo 
dt e iult J L (k,t), 
■oo 

/+oo 
dt e iuJt J T (k,t). (4) 
-oo 

In the simulation, Eqs. (T4j) are evaluated by using discrete Fourier transform in a period of 
655^"^, and the results are averaged over a longer period of 10480^. 
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III. RESULTS AND DISCUSSIONS 



In this section, we present our main simulation results for wave spectra, together with 



analytical results for dispersion relations for DAW 



25|, DTW [23|] , and those obtained by 



using QLCA [18( and HA [15|. The relevant theoretical derivations are briefly reviewed for 
the sake of completeness. 



DAW 



Let us begin with a weakly coupled state. We use here a fluid description for a mono- 
layer of dust particles levitating in plasma, without considering details of their mutual 
interactions. Assuming that a cold, 2D dust fluid occupies the plane z = in a Cartesian 
coordinate system with R = {x,y, z}, we let a d (r,t) and u d (r, t) be, respectively, number 
density per unit area and velocity field (having only the x and y components) of the dust 
fluid at the position r = {x, y} and at time t. The continuity equation and the momentum 



equation for the fluid are, respectively 
da d (r,t) 



25] 



dt 



+ V r Mr,t)u d (r,t)] = 0, (5) 



fM * ,UV/) + u d (r,t) • V||U d (r,t) = ^V||$(R,t) - vu d (r,t) (6) 



dt " in 



2=0 



where v is the Epstein drag coefficient. Note that the spatial differentiation in Eqs. ([5]) 

. d A d 

and (El) only includes tangential directions, viz., Vii = x- — h y^r- • The first term in the 

ox ay 

right-hand side of Eq. (jH]) indicates that, although the total electrostatic potential $(R, t) 
depends on all three spatial coordinates R = {r, z}, only the x and y components of the 
electrostatic force, evaluated in the plane z = 0, affect the motion of the dust fluid. The full 
spatial dependence of the electrostatic potential $ is determined by the Poisson equation in 
3D, 

V 2 $(R, t) = -4vre MR,t) - n e (R,t) - Z d a d (r,t) 5{z)} , (7) 
. d 

where V = Vii + z— . The electron and ion volume densities are given by Boltzmann 

oz 

relations, n e = n exp(e$/T e ) and rij = n exp(— e$/Tj), respectively, owing to the fact that 
the dynamics of massive dust particles is so slow that both electrons and ions are considered 
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to have enough time to reach their respective local equilibria, with uq being the equilibrium 
plasma density and T^ e \ the ion (electron) temperature (in energy units). 



The above equations can be solved perturbatively 
the 2D dust fluid in the form 



19 



25], giving a dielectric function of 



e{k, to) 

We define 

u 2 (k) = 



1 - 



2%e 2 Zja d0 



k 2 M 



1 



(8) 



2ire 2 Zjcr d0 (X D k) 2 
m\ D ^k 2 \ 2 D + 1 



m 



(9) 



and note that <p(k) is the 2D Fourier transform of the Yukawa potential. A dispersion 



relation for acoustic waves in the 2D dust fluid 



25j can be obtained from 



u(u + iu)=u%(k), (10) 

in analogy to the 3D case , |23j. The real part of the dispersion is shown in Fig. 

[5] (heavy dashed lines), while the corresponding discussion is postponed to the following 
subsection. 



B. Dust thermal wave (DTW) 



Note that, in the above derivation, we have assumed that the dust fluid is cold and, 
consequently, we neglected the direct thermal effect. Nevertheless, it would be interesting 
to see how this effect changes wave dispersion in the dust fluid. This effects can be easily 
retained by including the ideal-gas part of pressure P d k in the momentum balance equation, 
Eq. as follows [23 1 



fAMr /i + u d (r,t) ■ V||U d (r,t) = ^Vn$(R,t) 



/ \ 1 
- vu d (r,t) H 

z=o ma d0 



ViiP, 



(Ik ■ 



dt 1 in 

where P d k = r ) d k,BT d a d is the kinetic part of the pressure in dust component, and 7^ = 2 is 



the adiabatic index for the 2D dust system 



23] 



With this new momentum balance equation, the dielectric function of the 2D dust fluid 
becomes 



1 - 



ujiuj + iu) 



(12) 




FIG. 2: (Color online) (Color bar in arbitrary unit) Wave spectra of longitudinal mode for k = 1, 
v = O.Olcjpd and different T (nonideal gas states). Dashed lines are dispersion relations of DAW and 
result of semi-analytic approximation, respectively, while dash-dotted lines are results of DTW. 

with v t h = \j2kBTd/m being the thermal speed of dust particles, which gives the dispersion 
relation for DTW 

2 

uo{uo + iv)=uol{k) + ^^k 2 . (13) 
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The real part of the DTW dispersion is shown in Fig. [2] (dash-dotted lines) . 

One expects that the DAW and DTW modes, derived above, should be dominant in 
weakly coupled systems. However, to the best of our knowledge, this conjecture has not 
been examined in previous simulations. To elucidate this issue, we show in Fig. [2] wave 
spectra from simulations, together with the dispersion relations for both DAW and DTW 
for very low coupling strengths, at which systems are often regarded as being in a non- 
ideal gaseous state. One sees that, at such high temperatures, the collective modes are 
heavily damped at short wavelengths. There are essentially no collective modes beyond 
ka = 2 for T = 1, 2 and 5 and, moreover, these modes diminish at higher temperatures, 
or lower Ts. One would attribute this to Landau damping and viscous/collisional damping. 
Comparison with the above analytical results at T — 1 shows that the agreement between 
the DTW and the simulation is remarkably good, whereas the DAW displays noticeable 
discrepancy with simulation, indicating that the thermal effect is significant. With increasing 
T or, equivalently, decreasing temperature, a discrepancy between the DTW and simulation 
develops and becomes noticeable at, e.g., T = 5, whereas the DAW seems to agree better 
with simulation. As we'll see in the next subsection, this seemingly better agreement between 
DAW and simulation at T = 5 is just a coincidence, because the direct thermal effect and 
the strong coupling effect cancel each other, as will be discussed in the following subsection. 
The discrepancy arises immediately when T further increases, say at T = 10. 



C. Semi-analytic approximation 



The DAW and DTW derived in the previous subsections are essentially based on the 



mean-field theory of RPA 



22 



231 ]. in which the short-range interactions between dust par- 



ticles are neglected. Therefore, those results are applicable, in principle, only in weakly 
coupled situations, i.e., when T << 1. For strongly coupled systems, the short-range inter- 
particle interaction becomes important and it is necessary to take this effect into account. 
Following the line of reasoning in our fluid description above, this can be done immediately 
w introducing the interaction part of the pressure P d i into the momentum balance equation 



25|,|47| 



^ (r ' t) + u d (r,t).V l| u,(r,t) = ^V ll $(R,t) 



dt 



rn 



-vu d (r,t)-\ —V\\P dk ^ 

z =o ma d0 



1 



ma d0 



-ViiP*, (14) 
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where P d i usually contains information of system structure and interaction. A semi-analytic 



15 



results of dielectric function 
energy e c through density functional theory [47]. i. <\. 

'SP* 



25j can then be obtained by relating P d i to correlation 



VuP di 



5a d 



V\\a d . 



(15) 



Here one sees that a = (5Pu/Sa d )T is actually the isothermal compressibility [29|, |47j] and 
can be further written as 25 1 

crdo d 2 



a 



m do 2 ,. 



■[cr do e c (o do )]. 



(16) 



(See Ref. 25j and references therein for more details of the derivation.) Empirical 



of e r for 2D Yukawa system is now available through computer simulation 



15 



expression 



25] 



481 ] . One has 



e (k, uj) 



1 - 



lu(lu + iv) — ak 2 
and consequently the dispersion relation 



0;2( fc ) + mUtf 



(17) 



uoyuo + iv) 



co 2 (k) 



ak . 



(18) 



Here the term with a in Eq. (11711 can be regarded as a local field correction in the mean- 

[ i 

field theory accounting for the correction due to the strong coupling effect [47J. Because a 



is related to the correlation energy e c , which in turn can be obtained from simulation [48], 
this dispersion relation was entitled as a semi-analytic approximation [l5]. Typical value of 
a is about 0.2 ~ 0.3uj pd a. 

It should be noted that in Eq. (1181) . the direct thermal effect and the correlation effect 
have opposite influence on the dispersion relation. At very high temperature, say T = 1 
in Fig. [21 the correlation effect is negligibly small, so we saw a perfect agreement between 
DTW and simulation. With the increase of T, the former becomes smaller and smaller 
while the latter becomes more and more pronounced and at T = 5 the two effects somehow 
cancel each other, so we saw a good agreement between DAW and simulation. As can be 
expected, with the further increase of T, the correlation effect becomes dominant, so that 
the semi-analytic result gives a good agreement with simulation at F = 10 in Fig. [3 The 
similar effect had also been observed in GH theory of Kaw and Sen 29| for 3D dusty plasma 
liquids. 
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D. QLCA 




FIG. 3: (Color online) (Color bar in arbitrary unit) Wave spectra of Longitudinal mode for k = 1, 
r = 100 and v = 0.01u p d. Dashed lines are dispersion relations of DAW and QLCA and semi- 
analytic approximation. 

In previous subsection, we had come up to the concept of local field correction (LFC) and 
introduced the LFC through density functionals of interaction part of pressure, i. e., Eq. 



T6l) . In the following, we adopt the LFC in terms of the QLCA [18|, |26|, |35j, which had been 



"^ff-' ' n """' 4 " , m,, " OT ; "" , ""' ,,v " *' ."«■■■'->■,. i. * 



18l . Il9l . |35|. The corresponding dielectric function (for both longitudinal and transverse 



modes) of a strongly coupled 2D dusty plasma reads as follows 

SL/T (K U) = 1 - ^,"f (fc) n (19) 

uj{oj + iv) - D L / T (k) 

where Di{k) and Dx{k) are projections of the dynamical matrix in the longitudinal and 
transverse directions, respectively, and are functionals of the equilibrium radial distribution 
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FIG. 4: (Color online) (Color bar in arbitrary unit) Wave spectra of transverse mode for T = 100 
and 200, k = 1 and u = 0.01w p rf. Dashed lines are dispersion relation of QLCA. 



function, g(r). Detailed expressions for Di(k) and Df{k) can be found in Refs. 18|, l39 |. 
Note here that the direct thermal effect is neglected in QLCA and will be resumed in next 
subsection. 

The longitudinal and transverse modes are determined from these dielectric functions by 



14 



letting, 



e L {k,u) = 0, and e^ 1 {k,uj) = 0, (20) 
which give 

w(w + iu) = u 2 (k) + D L (k) (21) 

u(uj + iu) = D T (k), (22) 

for the longitudinal and transverse dispersion relations, respectively. 

Figure [3] shows a longitudinal phonon spectrum for T = 100 and v = 0.01u p d. Also 
shown are the dispersion relations of DAW, semi-analytic approximation and QLCA. A 
good agreement between the simulation and the QLCA is clearly seen in a broad range of 
wavelength on the long wavelength side, as expected [la, H9[. Discrepancy appears beyond 
the first Brillouin zone, especially for very large wave numbers, where the simulation displays 
a raising tail due to the DTW mode, while the QLCA converges to the Einstein frequency 
[lsI . Q] . This discrepancy is a consequence of the neglect of the direct thermal effects in the 
QLCA, which prompts us to discuss a suitable amendment to this theory, to be described 
in the following subsection. It is not surprising to find that the result of semi-analytic 
approximation is essentially the same as QLCA in the l ong wavelength region (up to around 
ka = 2.0), because actually a = \im k ^ D L (k)/k 2 [18|, |26|, |35| . However, the semi-analytic 
approximation breaks down for short wavelength. 

Figure H] shows the wave spectra of a transverse mode with v = 0.01u p d, for T = 100 
and 200, along with a dispersion relation from the QLCA. The transverse mode is heavily 
damped, especially at short wavelengths and in the system at a higher temperature, as can 
be seen from Fig. HI For T = 100, there are practically no transverse modes beyond ka ~ 6 
(implying a short-wavelength cutoff) while, for T = 200, the spectrum becomes very noisy 
beyond ka m 6, so that it is difficult to distinguish peaks due to collective motion. One 
can also notice that, for F = 100, the peak contour of the wave spectra does not seem 
to reach ka = at ui = 0, indicating a long- wavelength cutoff [311 ] . (This feature will be 
further discussed below.) Agreement between the simulation and the QLCA does not appear 
to be very satisfactory, neither at short nor at long wavelengths. The agreement at long 
wavelengths is somewhat improved for T = 200, in which case the long wavelength cutoff, 
implied by the simulation, vanishes. It should be noted that the melting point for this 

15 



system occurs at T* « 180 for k = 1 [18]. So, at T = 200, particles are almost frozen and 
wave propagation becomes anisotropic at short wave lengths. The task of resolving angular 
dependence of the wave dispersion lies beyond the capability of QLCA, but this issue can 
be tackled by the HA, as shown later on. 



E. QLCA with DTW in longitudinal mode 



2.0 



r=io../^pp^ Lmode 




FIG. 5: (Color online) (Color bar in arbitrary unit) Wave spectra of longitudinal mode for k = 1 
and v = O.Olujpd but different T covering both liquid and solid states. Thin dashed lines are results 
of DAW, while the heavy dashed lines are extensions of QLCA with DTW modes, i. e., EQLCA 



It is known 



la . l26l . |35| that the QLCA neglects the direct thermal effect, which is 
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responsible for the actual motions and migration of dust particles, and which gives rise to 
the so-called Bohm-Gross term (oc k 2 v 2 h ) in the third-frequency-momentum sum rule, and 



consequently in the longitudinal dispersion relation |26|]. Since this term is of the order of 
0(r _1 ) in comparison with terms due to the correlation effect, it is a good approximation to 
neglect it for T » 1 and, especially, at long wavelengths, ka « 1. At short wavelengths, 
ka >> 1, the k 2 factor could elevate the direct thermal effect up to a significant level. 
Nevertheless, this effect was not examined in previous simulations, and is therefore of interest 
for the present study. 

A simple extension of the QLCA, which would include the effect of direct thermal motion, 
can be made at the phenomenological level by taking into account a contribution from the 
kinetic part of pressure. Mathematically, this can be done by inserting the QLCA-LFC 
into Eq. ffl3|) . in which case it gives a correction to the longitudinal mode from the QLCA 
dielectric function, as follows 

1 



e (k,tu) = 1 — 



(23) 



lu(lu + w) — Di(k) 
and consequently the longitudinal dispersion relation becomes 

2 

u(u + iv) = cu 2 (k) + D L (k) + ^k 2 . (24) 

Figure [5] shows wave spectra of the longitudinal mode for different values of T covering 
a wide range in both liquid and solid states. Dispersion relations of both the DAW and the 
QLCA with the above DTW extension (denoted as EQLCA), i. e., EqfSU are also shown 
for comparison. It is seen that the raising DTW tails in the spectra become less and less 
significant with increasing T. The DTW mode is still noticeable even at T — 800, although 
it is very weak there, as shown. On the other hand, one also sees that the EQLCA captures 
the DTW tails in the spectra, in addition to exhibiting a consistently good agreement 
with simulation at small wave-numbers. Nevertheless, discrepancies still exist in some fine 
structures at intermediate wavelengths. Some of those, occurring for T > T*, are due to the 
angular-dependence of wave propagation in the solid state. 

It should be noted here that our extension is essentially based on mean-field theory of 
Rao 



23j and that QLCA is used only as an LFC. Therefore from this point of view, our 
extension is not strictly self-consistent and might not satisfy the third frequency sum rule. 
A more self-consistent theory, which can satisfy the third frequency sum rule, is now under 
developing and will be presented elsewhere 



49]. 
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F. QLCA with critical wave cutoff in transverse mode 



1000 



2000 



3000 




FIG. 6: (Color online) (Color bar in arbitrary unit) Transverse mode for T = 100, k = 1 and 
v = O.OliOpd- Dashed line is the fit according to Eq. (|27p . i. e., the extension of QLCA with cutoff 
wave- number. This gives a j/dm w 0.23uj p d and a cutoff wave-number of k c a = 0.48, as is labeled 
in the figure. The inserted figure shows k c a for k = 1, v = O.OlWp^ but different T. 



An another weakness of the QLCA lies in its inability to account for the diffusional and 
other damping effects that preclude the existence of long wavelength transverse waves in the 
liquid state fl8 |. Kalman et al. [3] provided a work-around by introducing a phenomeno- 
logical damping vdm- This gives rise to a small change in the transverse dielectric function 
of the QLCA, 



e (k,uj) 



oj{oj + iv + ivou) - D T {k) 



(25) 
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and consequently in the transverse dispersion relation, 



u(u + iv + wdm) — D T (k). (26) 
By assuming v = + , the real part of the dispersion relation can be written as 



Wr = \l Dr(k) — -^p-, (27) 

so that the condition Dx{k) — VdmK^) — gives rise to a critical cutoff wave-number k c . 

Note that vdm, and consequently k c , can be determined by directly fitting the simulation 
results to the above form [18|. Results of such fitting are shown in Fig. 0, along with a 
magnification of Fig. H] at small ks, where one can clearly notice the existence of a long 
wavelength cutoff k c . As a consequence, one can notice an improved agreement between the 
simulation and the QLCA with the extension introducing the cutoff. This extension of the 
QLCA does not affect much the dispersion at large k, where it behaves in much the same 
way as in Fig. H] because Udm ~ 0.23uj p d is relatively small. From the figure's inset, one 
sees that the cutoff wave-number k c increases with decreasing T. This tendency agrees with 



previous experimental observations 121 ]. 



G. HA for DLW 

In a perfect crystalline state, all particles are located at the points of a triangular lattice 
and they perform thermal oscillations around their equilibrium positions. In this case, one 
needs to adopt an another strategy to obtain dispersion relation for the DLW. This dispersion 
relation can be determined from the following eigenvalue problem [ijj] : 

||cu 2 (k)I-M(k)|| = 0, (28) 

where I is the 2D unit matrix. (Note that, although we have neglected the phenomenological 
damping in Eq. (1281) . it can be easily retained by replacing u 2 with uj{uj + 27).) The matrix 



M(k) is the interaction matrix, given by 15] 



i 

where a and (3 = x, y, and the summation over i includes all points on the triangular lattice. 
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FIG. 7: (Color online) (Color bar in arbitrary unit) Wave spectra of both longitudinal and trans- 
verse modes for n = 1, V = 1000 and v = 0.01uj p d- The upper panel shows the result of BD versus 
EQLCA and also DAW, while the lower panel for BD versus HA. The results of HA cover a whole 
period of ir/6 with an interval of 3 degree. 



For the HA 15], results are quite similar to those of QLCA. There are two branches given 



by, 



UJ 



L/T 



(M) 



[M xx + M x 



mil 



± - x l[M xx + M yy f + A{Ml y -M xx M yy) . 



(30) 



where M a p is given by Eq. (12U|) . and 6 is the polarization angle [151 ] . The subscripts L and T 
denote the longitudinal and transverse modes, and they correspond to + and — signs in Eq. 
f |30|) . respectively. Therefore, wave propagation depends on the angle 9, and its dispersion 
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relations are functions with the period of 7r/6 due to the hexagonal symmetry [15j. The 
angular dependence arises because the system is now anisotropic at short wavelengths. 

Figure [7| shows wave spectra of both longitudinal and transverse modes in a crystalline 
state with T = 1000. Also shown are the dispersion relations of the HA, EQLCA and DAW. 
The HA curves cover the whole period of 7r/6 with an increment of 3 degrees. It is seen 
that the simulation agrees very well with the HA, and that the polarization effects in the 
simulation spectra are fully captured by the HA. Good agreement between the EQLCA and 
the simulation is retained in the first Brillouin zone. For large wave-numbers, the EQLCA 
cannot resolve the polarization, i.e., the angular-dependence of the dispersion. Nevertheless, 
it captures the correct tendency for oscillations. Moreover, at this range of T values, the 
QLCA regains its credibility because both the thermal effect and the long wavelength cutoff 
corrections become unimportant, and the EQLCA essentially reduces to the standard QLCA. 
At the same time, it is remarkable to see that the mean field DAW is so robust, even in 
a crystalline state, where it still offers a fairly good agreement with the simulation at long 
wavelengths. 



H. Damping effect 

We finally investigate the effects of damping on wave spectra. Generally speaking, there 
are two types of damping mechanisms for collective modes, i.e., the intrinsic and external 
ones. The former type includes dampin g du e to diffusion, viscosity, thermal conduction 



and, possibly, also Landau-like damping 



26. 



46 



471 ]. while the latter type mainly comes 



from the neutral gas. The external damping is simply given by the imaginary part of the 
dispersion relation by u>i = —v/2, , whereas the intrinsic damping of collective modes involves 
complicated physical processes, and so far there has been no well founded theory of such 
processes for SCDPs. Development of such a theory may be an interesting topic for future 
research, while, at this stage, we present several results from our simulation. 

Figure [8] shows the profiles, or, more precisely, the cross sections at several fixed wave 
numbers, of the longitudinal wave spectra, £(k, u), for T = 100 and with different neutral 
damping rates v. These profiles exhibit typical resonance shapes which describe the spectral 
distributions shown in Figs. 1-7, while their peak positions in the u-k plane correspond to the 
dispersion relations. The above mentioned analytical theories give only the peak positions 
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pd 



FIG. 8: Profiles of longitudinal wave spectra for certain wavenumbers with k = 1, T = 100 but 
different damping rate v. 



in the u-k plane. In reality, these peaks are never delta-like functions, but rather have 
finite amplitudes and consequently finite widths due to various damping effects, as is shown 
in Fig. [HI Physically, the peak amplitudes indicate the amount of energy concentrated in 
various modes, while the widths of the peaks at half heights (denoted as the half-widths, 
Aco>i/2) are related to damping 9]. One notices a general tendency that the peak amplitudes 
decrease with the increase of wave-number. Similar tendency was also observed in classical 
MD simulation [35]. In addition, with the increase of the neutral gas damping rate u, all 
the heights decrease, while the widths increase. These tendencies are depicted more clearly 
in Fig. [91 which shows the amplitudes and the half-widths of C(k,u) and T(k, u) versus 
wave-number for different v. 
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FIG. 9: Peak-amplitudes and half-widths of both longitudinal and transverse wave spectra 

for k = 1, T = 100 but different neutral gas damping rate v. 

IV. CONCLUSION 

By using the Brownian dynamics simulation, we have investigated the wave spectra of 
2D dusty plasmas in states that cover a full range from a non-ideal gas to crystalline. The 
results are critically compared with dispersion relations for the dust-acoustic wave, dust- 
thermal wave, and those from the quasi-localized charge approximation and the harmonic 
approximation. In particular, simply extensions are considered to include the direct thermal 
effect and the critical cutoff wavenumber into the QLCA. It is found that, for a non-ideal 
gaseous state (e.g., with T = 1), the longitudinal DTW is in a remarkably good agreement 
with the simulation, while for a perfect crystalline state (e.g., T = 1000), the HA agrees 
very well with the simulation. In the states between those two opposing ends, an overall 
good agreement between the extended versions of the QLCA and the simulation was found 
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for a wide range of wavelengths. The damping effect is also briefly discussed in terms of the 
peak amplitudes and the half-widths of the current-current correlation function for different 
neutral gas damping rate v. A general tendency is that the peak amplitudes decrease with 
the increase of wave-number and the increase of u, while the widths increase during these 
courses. 
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